// $Id: cylrrs.cpp 14 2008-02-23 15:12:07Z msaushkin $
/**
* @file
* Реализация пространства имён cylrrs, используемого для восстановления полей остаточных напряжений и 
* деформаций в цилиндрическом образце.

* Компиляция модуля осуществляется командой (в GNU/Linux):
* g++ -Wall -I/usr/local/include -I. -c cylrrs.cpp
* @author Саушкин М.Н.
*/

#include<iostream> 
#include "cylrrs.hpp"
#include<math.h> 
#include<gsl/gsl_sf_erf.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_pow_int.h>

/**
* Структура, описывающая основные параметры математической модели класса, которые задаются
* открытыми атрибутами. По задаваемым атрибутам (параметрам) определяются остальные.
* Структура используется  в пространстве имён cylrrs.
*/
struct cyl_params{
                  double radius;  //!< радиус цилиндрического образца, мм (задаваемый параметр).
                  double h0;      //!< глубина обнуления эпюры \f$\sigma_\theta^{\rm res}\f$, мм (задаваемый параметр).
                  double sigmazv; //!< значение эпюры \f$\sigma_\theta^{\rm res}\f$ на поверхности \f$r=a\f$, кг/мм\f$^2\f$  (задаваемый параметр).
                  double yuong;  //!< модуль Юнга материала из которого изготовлен образец, кг/мм\f$^2\f$  (задаваемый параметр).
                  double nu;     //!< коэффициент Пуассона \f$\nu\f$, (задаваемый параметр).
                  double alpha;  //!< "параметр упрочнения" или "гипотеза упрочнения" (задаваемый параметр).
                  double b;      //!< параметр аппроксимации эпюры \f$\sigma_\theta^{\rm res}\f$, мм (вычисляемый параметр).
                  double sigma0; //!< параметр аппроксимации эпюры \f$\sigma_\theta^{\rm res}\f$, кг/мм\f$^2\f$  (вычисляемый параметр).
                  double sigma1; //!< параметр аппроксимации эпюры \f$\sigma_\theta^{\rm res}\f$, кг/мм\f$^2\f$  (вычисляемый параметр).
                  double epsz0;  //!< величина полной осевой деформации цилиндра (вычисляемый параметр).
               };

/**
* Пространство имён для описания и реализации математической модели восстановления полей остаточных 
* напряжений и пластических деформаций в цилиндрическом образце.
*/
namespace cylrrs{

   /** Основные параметры математической модели класса.
   */
   cyl_params params;
 
   /** Инициализация параметров модели.
   */
   int init(double params_radius, double params_h0, double params_sigmazv, double params_yuong, double params_nu, double params_alpha);

   /**
   * Функция для нахождения параметра модели \f$b\f$.
   *
   * Этой функцией задаётся уравнение для нахождения параметра \f$b\f$, которое решается методом 
   * дихотомии. 
   *
   * @param b
   *    Параметр аппроксимации эпюры \f$\sigma_\theta^{\rm res}\f$, мм 
   * @param radius
   *    Радиус цилиндрического образца, мм
   * @param h0
   *    Глубина обнуления эпюры \f$\sigma_\theta^{\rm res}\f$, мм
   * @return 
   *   Значение функции при заданных параметрах. Если функция равна 0, то параметр \f$b\f$ искомый.
   */
   double f_cyl(double b, double radius, double h0);

   /** Вычисление параметра \f$b\f$ аппроксимации \f$\sigma_\theta^{\rm res}\f$.
   */
   double b_Calc(double (*fun) (double, double, double), double, double, double, void *);

   /** Вычисление параметра \f$\sigma_1\f$ аппроксимации \f$\sigma_\theta^{\rm res}\f$.
   */
   double sigma1_Calc(void * );

   /** Вычисление параметра \f$\sigma_0\f$ аппроксимации \f$\sigma_\theta^{\rm res}\f$.
   */
   double sigma0_Calc(void * );

   /** Вычисление \f$\sigma_\theta^{\rm res}(r)\f$ по радиусу.
   */
   double sigmatheta(double, void * );

   /** Вычисление \f$\sigma_r^{\rm res}(r)\f$ по радиусу.
   */
   double sigmar(double, void * );

   /** Подынтыгральная функция в выражении для \f$q_\theta(r)\f$.
   */
   double qtheta_integrand(double, void *);

   /** Вычисление \f$q_\theta(r)\f$ по радиусу.
   */
   double qtheta(double, void * );

   /** Вычисление \f$q_z(r)\f$ по радиусу.
   */
   double qz(double, void * );

   /** Вычисление \f$q_r(r)\f$ по радиусу.
   */
   double qr(double, void * );

   /** Подынтыгральная функция в выражении для \f$\varepsilon_z^0\f$.
   */
   double epsz0_integrand(double, void *);

   /** Вычисление параметра \f$\varepsilon_z^0\f$.
   */
   double epsz0_Calc (void * );

   /** Вычисление \f$\sigma_z^{\rm res}(r)\f$ по радиусу.
   */
   double sigmaz(double r, void * p);

   double getsigmaz(double);
   double getsigmar(double);
   double getsigmatheta(double);
   double getqz(double);
   double getqr(double);
   double getqtheta(double);
   double getradius();
   double geth0();
   double getsigmazv();
   double getyuong();
   double getnu();
   double getalpha();
   double getb();
   double getsigma0();
   double getsigma1();
   double getepsz0();
   void paraminfo();
};


double cylrrs::f_cyl(double b, double radius, double h0){
   return exp(-gsl_pow_2(h0/b))-b*sqrt(M_PI)*gsl_sf_erf(radius/b)/(2.0*radius); 
   };

int cylrrs::init(double params_radius, double params_h0, double params_sigmazv, double params_yuong, double params_nu, double params_alpha){
   params.radius = params_radius;
   params.h0 = params_h0;
   params.sigmazv = params_sigmazv;
   params.yuong = params_yuong;
   params.nu = params_nu;
   params.alpha = params_alpha;
   params.b = b_Calc(f_cyl, .00000001, .3, 1e-12, &params);
   params.sigma1 = sigma1_Calc(&params);
   params.sigma0 = sigma0_Calc(&params);
   params.epsz0  = epsz0_Calc(&params);
   return 1;
};

double cylrrs::b_Calc(double (*fun) (double, double, double), double a, double b, double epsilon, void * p){
  struct cyl_params * here_params = (struct cyl_params *)p;
  double radius = (here_params->radius);
  double h0    = (here_params->h0);
  double c;
  while (b-a > epsilon){
     c = (a+b)/2.0;
     if(fun(b, radius, h0) * fun(c, radius, h0) < 0)
       a = c;
     else
       b = c;
    }
  return (a+b)/2.0;
};

double cylrrs::sigma1_Calc(void * p){
  struct cyl_params * here_params = (struct cyl_params *)p;
  double h0      = (here_params->h0);
  double b       = (here_params->b);
  double sigmazv = (here_params->sigmazv);
  return sigmazv/(exp(-gsl_pow_2(h0/b))-1.0);
};

double cylrrs::sigma0_Calc(void * p){
  struct cyl_params * here_params = (struct cyl_params *)p;
  double sigmazv = (here_params->sigmazv);
  double sigma1  = (here_params->sigma1);
  return sigmazv + sigma1;
};

double cylrrs::sigmatheta(double r, void * p){ 
  struct cyl_params * here_params = (struct cyl_params *)p;
  double radius = (here_params->radius);
  double sigma0 = (here_params->sigma0);
  double sigma1 = (here_params->sigma1);
  double b      = (here_params->b);
  return sigma0 - sigma1*exp( -gsl_pow_2(radius-r) /gsl_pow_2(b) );
};

double cylrrs::sigmar(double r, void * p){
  struct cyl_params * here_params = (struct cyl_params *)p;
  double radius = (here_params->radius);
  double sigma0 = (here_params->sigma0);
  double sigma1 = (here_params->sigma1);
  double b      = (here_params->b);
  if (r!=0) 
    return sigma0 - (sigma1*b*sqrt(M_PI))*(gsl_sf_erf(radius/b)-gsl_sf_erf((radius-r)/b))/(2.0*r);
  else
    return sigma0;
};

double cylrrs::qtheta_integrand(double z, void * p){
  struct cyl_params * here_params = (struct cyl_params *)p;
  double alpha = (here_params->alpha);
  double nu    = (here_params->nu);
  double l = 1+alpha-alpha*nu;
  double k = 1+alpha*nu;
  double m = l/k;
  return pow(z, m)*(sigmar(z, here_params)+(1+alpha)*sigmatheta(z, here_params) );
};

double cylrrs::qtheta(double r, void * p){
   struct cyl_params * here_params = (struct cyl_params *)p;
   double alpha = (here_params->alpha);
   double nu    = (here_params->nu);
   double yuong = (here_params->yuong);
   double a =  (1+nu)*(1-2*nu)/(yuong*(1+alpha*nu)*(1+alpha*nu));
   double b =  (1+nu)/(yuong*(1+alpha*nu));
   double c = -(2+alpha)/(1+alpha*nu);

   if (r!=0){
         gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);
         double IntRes, result, error;
         gsl_function F;
         F.function = &qtheta_integrand;
         F.params = here_params;
         gsl_integration_qags(&F, 0, r, 0, 1e-10, 1000, w, &result, &error);
         IntRes = result;
         gsl_integration_workspace_free(w);
         result = pow(r, c)*a*IntRes - b*((1-nu)*sigmatheta(r, here_params)-nu*sigmar(r, here_params));
         if (fabs(result)<1e-19) 
            return 1e-20;
         else return result;
         }
      else return 0.0;
};

double cylrrs::qz(double r, void * p){
   struct cyl_params * here_params = (struct cyl_params *)p;
   double alpha = (here_params->alpha);
   return alpha*qtheta(r, here_params);
};

double cylrrs::qr(double r, void * p){
   struct cyl_params * here_params = (struct cyl_params *)p;
   double alpha = (here_params->alpha);
   return -(1 + alpha)*qtheta(r, here_params);
};

double cylrrs::epsz0_integrand(double r, void * p){
   struct cyl_params * here_params = (struct cyl_params *)p;
   double nu    = (here_params->nu);
   double yuong = (here_params->yuong);
   double a = nu/yuong;
   return r*( qz(r, here_params)-a*( sigmar(r, here_params) +  sigmatheta(r, here_params)));
};

double cylrrs::sigmaz(double r, void * p){
   struct cyl_params * here_params = (struct cyl_params *)p;
   double epsz0 = (here_params->epsz0);
   double yuong = (here_params->yuong);
   double nu    = (here_params->nu);
   return yuong*(epsz0-qz(r, here_params))+nu*(sigmar(r, here_params) +  sigmatheta(r, here_params));
};


double cylrrs::epsz0_Calc(void * p){
   struct cyl_params * here_params = (struct cyl_params *)p;
   double radius  = (here_params->radius);
   double IntRes, result, error;
   gsl_integration_workspace * w = gsl_integration_workspace_alloc(1000);
   gsl_function F;
   F.function = &epsz0_integrand;
   F.params = here_params;
   gsl_integration_qags(&F, 0, radius, 0, 1e-10, 1000, w, &result, &error);
   IntRes = result;
   gsl_integration_workspace_free(w);
   return 2*IntRes/(radius*radius);
};

double cylrrs::getsigmaz(double r){
   return sigmaz(r, &params);
};

double cylrrs::getsigmar(double r){
   return sigmar(r, &params);
};

double cylrrs::getsigmatheta(double r){
   return sigmatheta(r, &params);
};

double cylrrs::getqz(double r){
   return qz(r, &params);
};

double cylrrs::getqr(double r){
   return qr(r, &params);
};

double cylrrs::getqtheta(double r){
   return qtheta(r, &params);
};

void cylrrs::paraminfo(){
   std::cout << "radius = " << params.radius << " mm" << "\n";
   std::cout << "h0     = " << params.h0 << " mm" << "\n";
   std::cout << "sigmazv= " << params.sigmazv << " kg/mm2" << "\n";
   std::cout << "yuong  = " << params.yuong << " kg/mm2" << "\n";
   std::cout << "nu     = " << params.nu << "\n";
   std::cout << "alpha  = " << params.alpha << "\n";
   std::cout << "b      = " << params.b << " mm" << "\n";
   std::cout << "sigma0 = " << params.sigma0 << " kg/mm2" << "\n";
   std::cout << "sigma1 = " << params.sigma1 << " kg/mm2" << "\n";
   std::cout << "epsz0  = " << params.epsz0 << "\n";
};

double cylrrs::getradius(){
   return params.radius;
}

double cylrrs::geth0(){
   return params.h0;
};

double cylrrs::getsigmazv(){
   return params.sigmazv;
};

double cylrrs::getyuong(){
   return params.yuong;
};

double cylrrs::getnu(){
   return params.nu;
};

double cylrrs::getalpha(){
   return params.alpha;
};

double cylrrs::getb(){
   return params.b;
};

double cylrrs::getsigma0(){
   return params.sigma0;
};

double cylrrs::getsigma1(){
   return params.sigma1;
};

double cylrrs::getepsz0(){
   return params.epsz0;
};
